﻿#include "qpcrModel.h"

using namespace QPCRModel;



#pragma region CALCULATION{


int Model::calculate(){

    _slopeList.clear();

    _Rmin = _rawPointList[3].y();
    _Rmax = 0;
    for(int c = _clipStart; c <= _cycleNum; c ++){
        _Rmax = std::max(_rawPointList[c].y(), _Rmax);
        _Rmin = std::min(_rawPointList[c].y(), _Rmin);
    }


    QVector<QPointF> trainList;
    if(_fixedBaseEnd < 0){          // 自动判断基线结束点
        // 1. 确定指数期位置 ================================================================================
        //_CSPE = calBaseline();
        //_CSPE = findExpStage();
        _CSPE = findExpStage2();

        if (_CSPE < _baseLen){
            if(_CSPE != -1){
                _resState = "===AnalysisResult: Failed- CSPE<baseLen";
            }
            return -1;
        }

        // 2. 参与拟合的点的点、初始化参数 =====================================================================================
        int train_start = _CSPE - _baseLen;
        int train_num = _cycleNum - train_start + 1;

        // 设置补偿点来调整基线期位置
        // 求baselen的平均值
        double baseAvg = 0;
        for(int i = train_start; i < train_start + _baseLen; i ++){
            baseAvg += _rawPointList[i].y();
        }baseAvg /= _baseLen;
        QVector<QPointF> compensationList;
        int compensation_num = 3;
        for(int i = compensation_num; i >= 1; i --){
            compensationList.append(QPointF(train_start - i, baseAvg));
        }
        trainList.append(compensationList);
        trainList.append(_rawPointList.mid(train_start, train_num));


        func->init(_CSPE, _CEPE, _cycleNum, _Rmin, _Rmax);    // 初始化
    }else{                          // 固定基线结束点
        if(_fixedBaseEnd >= _cycleNum){
            _resState = "===AnalysisResult: Failed- FixedBaseEnd < CycleNum";
            return -1;
        }
        _CSPE = _fixedBaseEnd;
        double baselineY = _rawPointList[_fixedBaseEnd].y();
        QVector<QPointF> compensationList;
        int compensation_num = 8;
        for(int i = compensation_num; i >= 1; i --){
            compensationList.append(QPointF(_fixedBaseEnd-i, baselineY));
        }
        trainList.append(compensationList);
        trainList.append(_rawPointList.mid(_fixedBaseEnd, _cycleNum));
        func->init(_CSPE, _CSPE+5, _cycleNum, baselineY, _Rmax);    // 初始化
    }


    double baseComp = 2000;         // 在拟合过程中给y加了一个补偿，防止在基线等于0的时候LM算法挂掉……
    func->y0 += baseComp;
    for(int i = 0; i < trainList.size(); i ++){
        trainList[i].setY(trainList[i].y()+baseComp);
    }
    if(LevenbergMarquardt(trainList)==-1){                // 拟合
        return -1;
    }
    func->y0 -= baseComp;


    // 3. 参数更新 ===================================================================================
    if(updateParam() == -1){
        return -1;
    }

    return 0;

}

int Model::calBaseline(){
    int c1 = 4, c2 = 15;
    for( ; c2 < _cycleNum-1; c2 ++){
        double k, b;
        //qDebug("尝试拟合直线%d到%d", c1, c2);
        double stderror = linearFitting(_rawPointList.mid(c1, c2-c1+1), &k, &b);

        if(k*_cycleNum < -(_Rmax-_Rmin)/15){
            continue;
        }
        int flag = 1;
        for(int c = c2+1; c <= _cycleNum; c ++){
            //qDebug("distancePL: %f", distancePL(_rawPointList[c], k, b));
            //qDebug("\t检测点%d", c);

            //if(distancePL(_rawPointList[c], k, b) <= stderror){
            if(- ((_rawPointList[c].x()*k - _rawPointList[c].y() + b) / sqrt(k*k+1)) <= stderror){
                flag = 0;
                //qDebug("\t %d不符合要求 重新拟合", c);
                break;
            }
        }
        if(flag == 1){
            _base_k = k;
            _base_b = b;
            _baselinePointList.clear();
            for(int c = 0; c <= _cycleNum; c ++){
                _baselinePointList.append(QPointF(c, k*c+b));
            }
            _RNoise = stderror;
            return c2+1;
        }
    }
    return -1;
}


int Model::findExpStage(){              // 已弃用

    // 指数期判断多加一次滤波
    QVector<QPointF> filterPointList = this->filter(_rawPointList, 0);

    QVector<QPointF> expPointList;
    QVector<double> corrList;
    for(int c = 0; c < _cycleNum - _expLen; c ++){
        int c1 = c, c2 = c + _expLen;
        double r_min = filterPointList[c1].y();
        for(int i = c1; i <= c2; i ++){
            r_min = std::min(filterPointList[i].y(), r_min);
        }
        expPointList.clear();
        for(int i = c1; i <= c2; i ++){
            expPointList.append(QPointF(filterPointList[i].x(), log10(std::max(1.0, filterPointList[i].y()-r_min))));
        }
        corrList.append(correlationPointList(expPointList));
        //qDebug("Correlation[%d %d]: %f", c1, c2, corrList[corrList.size()-1]);
    }

    int flag = -1;
    double maxCorr = -1;
    int    maxCorr_posi = -1;
    for(int i = corrList.size()-2; i > 3; i --){

        if(flag > 0){
            if(maxCorr < corrList[i]){
                maxCorr = corrList[i];
                maxCorr_posi = i;
            }
            flag --;
        }
        if(flag == 0){
            break;
        }

        if(corrList[i] > 0.8 && corrList[i+1] > 0.7 && corrList[i-1] > 0.7){
            flag = 5;
            maxCorr = corrList[i];
            maxCorr_posi = i;
        }
    }

    //qDebug("maxCorrPosi: %d", maxCorr_posi);

    return maxCorr_posi > 3 ? maxCorr_posi : -1;

}




int Model::findExpStage2(){


    // y轴归一化后的斜率
    QVector<QPointF> pointList = this->filter(_rawPointList, 0);
    // 计算Rmin、Rmax
    double Rmin = pointList[4].y(), Rmax = pointList[4].y();
    for(int c = _clipStart; c <= _cycleNum; c ++){
        Rmin = std::min(pointList[c].y(), Rmin);
        Rmax = std::max(pointList[c].y(), Rmax);
    }
    for(int c = 0; c <= _cycleNum; c ++){
        pointList[c].setY((pointList[c].y()-Rmin) / (Rmax-Rmin) * _cycleNum);
    }

    QVector<double> slopeList;
    for(int c = 0; c <= _cycleNum - _expLen; c ++){
        double k, b;
        this->linearFitting(pointList.mid(c, _expLen), &k, &b);
        slopeList.append(k);
    }
    _slopeList = slopeList;

    int spl = -1, spe = -1; // 线性期位置 / 指数期位置
    // 查找最大斜率位置
    for(int c = slopeList.size()-1; c >= _clipStart; c --){
        if(slopeList[c] > _slopeThr && spl == -1){
            spl = c;
            while(c --){
                if(slopeList[c] < 2){
                    spe = c;
                    break;
                }if(slopeList[c] > slopeList[spl]){
                    spl = c;
                }
            }
            break;
        }
    }
    if(spl == -1){
        _resState = "===AnalysisResult: Failed- SPL=-1";
        return -1;
    }
    // 查找是否有明显下降的区域
    for(int c = _clipStart; c < slopeList.size(); c ++){
        if(slopeList[c] < -slopeList[spl]*_slopeDownThr){
            _resState = "===AnalysisResult: Failed- slopeDownThr exceed limit";
            return -1;
        }
    }
    // 查找扩增结束位置
    _CEPE = _cycleNum + 1;
    for(int c = spl; c <= _cycleNum-_expLen; c ++){
        if(slopeList[c] < 2){
            _CEPE = c + _expLen;
            break;
        }
    }

    _CSPL = spl;
    _CSPE = spe;

    return spe;
}


// 扩增效率计算方法
double Model::calEfficiency(){
    /*
    (1) (Fn+1)/(Fn)
    (2) log(ΔRn)的斜率
    (4) 指数模型
    (5) 扩增效率模型（例如FPK-PCR）
    */
    double E;

    // 算法1
    E = 1;
    double exp1 = _CSPE, exp2 = std::min(_CSDM, (double)_cycleNum);
    int expLen = 0;
    for(double c = exp1; c < exp2; c ++){
        E *= (func->f(c+1)-_bgValue) / (func->f(c)-_bgValue);
        expLen ++;
    }

    E = pow(10, log10(E) / expLen);

//    // 算法2
//    QVector<QPointF> logExpPoints;
//    for(int c = _CSPE; c <= _CSDM; c ++){
//        logExpPoints.append(QPointF(c, log(func->f(c)-_bgValue)));
//    }
//    double ek, eb;
//    E = linearFitting(logExpPoints, &ek, &eb);

    return E;



}


int Model::LevenbergMarquardt(const QVector<QPointF> &dataList){


    _LM_logs.clear();

    int T = _iterationNum;
    double initMiu = 20;
    double miu = initMiu, last_L = -1;

    int psize = func->getParamNum();
    _LM_logs.push_back("DataList长度:" + QString::number(dataList.size()));
    _LM_logs.push_back(QString("Init: a: %1, b:%2, x0:%3, y0: %4, d: %5, miu:%6").arg(
                                                           QString::number(func->a),
                                                           QString::number(func->b),
                                                           QString::number(func->x0),
                                                           QString::number(func->y0),
                                                           QString::number(func->d),
                                                           QString::number(miu)));

    for(int t = 0; t < T; t ++){
        miu = miu > initMiu ? initMiu : miu;

        Eigen::VectorXd J(psize);                   J.setZero();
        Eigen::MatrixXd H(psize,psize);             H.setZero();
        Eigen::VectorXd g(psize);                   g.setZero();
        double loss2 = 0, loss = 0;
        for(int c = 0; c < dataList.size(); c ++){
            double xi = dataList[c].x();
            double yi = dataList[c].y();
            J += func->jacobian(xi, yi);
            loss += func->loss(xi, yi);
            loss2 += func->loss2(xi, yi);

            H +=  J * J.transpose();
            g += -loss * J;
        }

        //if(t % 10 == 0)
        //qDebug("time %d Loss: %.3f, a: %.3f, b:%.3f, x0:%.3f, y0: %.6f, d: %.3f, miu:%.3f", t, loss2, func->a, func->b, func->x0, func->y0, func->d, );
        _LM_logs.push_back(QString("time %1 Loss: %2, a: %3, b:%4, x0:%5, y0: %6, d: %7, miu:%8").arg(
                               QString::number(t),
                               QString::number(loss2),
                               QString::number(func->a),
                               QString::number(func->b),
                               QString::number(func->x0),
                               QString::number(func->y0),
                               QString::number(func->d),
                               QString::number(miu)));

        Eigen::VectorXd delta(psize);
        int while_num = 100;
        while(while_num --){
            Eigen::MatrixXd D(psize, psize), I(psize, psize);
            I = Eigen::MatrixXd::Identity(psize, psize);
            D = J * J.transpose();
            D = (I.array() * D.array()).matrix();
            delta = (H + miu * D).inverse() * g;
            for(int tmp = 0; tmp < delta.size(); tmp ++){
                if(delta[tmp] != delta[tmp]){
                    _resState = "===AnalysisResult: Failed- LM-delta irreversible";
                    return -1;
                }
            }
            Eigen::VectorXd newX(psize); newX = func->getV() + delta;
            double newLoss = 0;
            for(int c = 0; c < dataList.size(); c ++){
                double xi = dataList[c].x();
                double yi = dataList[c].y();
                newLoss += func->loss(xi, yi, newX);
            }
            double dl1 = newLoss - loss;
            double dl2 = J.transpose() * delta;
            double rho = dl1 / dl2;

            // 范围判断
            if(rho < 0.25){
                miu *= 2;
                continue;
            }
            if(1){
                int flag_wrong_param = 0;
                for(int i = 0; i < psize; i ++){
                    if(func->getV()[i]+delta[i] < func->getRange(i).first || func->getV()[i]+delta[i] > func->getRange(i).second){
                        miu *= 4;
                        flag_wrong_param = 1;
                        break;
                    }
                }
                if (flag_wrong_param == 1){
                    continue;
                }
            }
            if (rho > 0.75){
                miu /= 3;
                break;
            }

            break;
        }

        Eigen::VectorXd newX(psize);
        newX = func->getV() + delta;
        func->setV(newX);


        if(last_L != -1 && last_L / loss2 < 1.0000001 && loss2 / last_L < 1.0000001){
            break;
        }last_L = loss2;

    }

    return 0;
}

int Model::updatePointList(){
    _fittedPointList.clear();
    _fittedPointList_delta.clear();
    _fittedPointList_log.clear();
    _fittedPointList_norm.clear();
    _fittedPointList_log_delta.clear();

    if (_S_SHAPED){
        for(int c = 0; c <= _cycleNum; c ++){
            // 指数拟合处理
            if(func->f(c) > _Rmax){
                _fittedPointList.append(QPointF(c, _Rmax));
            }else{
                _fittedPointList.append(            QPointF(c, func->f(c)));
            }


            _fittedPointList_delta.append(      QPointF(c, func->f(c)-_bgValue));
            _fittedPointList_log.append(        QPointF(c, log10(std::max(10.0, func->f(c)))));
            _fittedPointList_log_delta.append(  QPointF(c, log10(std::max(10.0, func->f(c)-_bgValue))));
            _fittedPointList_norm.append( QPointF(c, (func->f(c)-func->f(1))/(func->f(_cycleNum)-func->f(1))));
        }
        // 计算拟合结果可信度
        double fittingError = 0;
        for(int c = std::max(_clipStart, (int)_CSPE-_baseLen); c <= _cycleNum; c ++){
            fittingError += sqrt( pow(_fittedPointList[c].x()-_rawPointList[c].x(), 2) + pow(_fittedPointList[c].y()-_rawPointList[c].y(), 2) );
        }
        fittingError /= _cycleNum - std::max(_clipStart, (int)_CSPE-_baseLen);
        _fittingConf = (_Rmax-_Rmin - fittingError) / (_Rmax-_Rmin);
    }
//    // 拟合失败数据处理 - 版本1
//    if(!_S_SHAPED || _fittingConf < _minConf){
//        _fittingConf = -1;
//        _fittedPointList.clear();
//        _fittedPointList_delta.clear();
//        _fittedPointList_log.clear();
//        _fittedPointList_log_delta.clear();
//        _fittedPointList_norm.clear();

//        _fittedPointList.append(this->filter(_rawPointList));
//        for(int c = 0; c <= _cycleNum; c ++){
//            _fittedPointList_delta.append(      QPointF(c, _fittedPointList[c].y()-_Rmin));
//            _fittedPointList_log.append(        QPointF(c, log10(std::max(10.0, _fittedPointList[c].y()))));
//            _fittedPointList_log_delta.append(  QPointF(c, log10(std::max(10.0, _fittedPointList[c].y()-_Rmin))));
//            _fittedPointList_norm.append(       QPointF(c, (_fittedPointList[c].y()-_Rmin)/(_Rmax-_Rmin)));
//        }
//    }

    // 拟合失败数据处理 - 版本2
    if(!_S_SHAPED || (_fittingConf < _minConf && _fixedBaseEnd==-1)){
        _fittingConf = -1;
        _fittedPointList.clear();
        _fittedPointList_delta.clear();
        _fittedPointList_log.clear();
        _fittedPointList_log_delta.clear();
        _fittedPointList_norm.clear();

        // 求原始数据均值
        double RAvg = 0;
        for(int c = _clipStart; c <= _cycleNum; c ++){
            RAvg += _rawPointList[c].y();
        }RAvg /= _cycleNum - _clipStart + 1;
        for(int c = 0; c <= _cycleNum; c ++){
            _fittedPointList.append(QPointF(c, RAvg));
            _fittedPointList_delta.append(QPointF(c, 0));
            _fittedPointList_log.append(QPointF(c, log10(std::max(RAvg, 10.0))));
            _fittedPointList_log_delta.append(QPointF(c, 1));
            _fittedPointList_norm.append(QPointF(c, 0));
        }
    }

    return 0;
}


int Model::updateParam(){

    // 注：各种参数的计算还是基于拟合后的函数，而不是数据点

    _bgValue = func->y0;
    if(_functype == 0){
        _CFDM = func->x0;
    }else{
        _CFDM = func->x0 + func->b * log(func->d);
    }
    //_RNoise = (func->g(_CFDM) - func->y0) / 20;
    _CSPE = func->g( _bgValue + (func->f(_cycleNum)-_bgValue)*0.03 );

//    // 修改CSPE算法，改为曲率最大值
//    double curvity_max = -1, CSPE;
//    for(double c = 1; c <= _cycleNum; c += 0.1){
////        double curvity = std::fabs(func->curvity(c, _cycleNum));
//        double curvity = func->curvity(c, _cycleNum);
//        if(curvity > curvity_max){
//            curvity_max = curvity;
//            CSPE = c;
//        }
//        qDebug() << "c =" << c << "\tcurvity =" << curvity;
//    }_CSPE = CSPE;
//    for(double c = CSPE-0.6; c < CSPE + 0.6; c += 0.05){
////        double curvity = std::fabs(func->curvity(c, _cycleNum));
//        double curvity = func->curvity(c, _cycleNum);
//        if(curvity > curvity_max){
//            curvity_max = curvity;
//            _CSPE = c;
//        }
//    }
    // 还是不改了


    // 数据筛选
    if(_fixedBaseEnd == -1){
        if(std::isinf(_CSPE) || std::isnan(_CSPE)){
            _resState = "===AnalysisResult: Failed- CSPE=NAN";
            return -1;
        }if(_CSPE < 2 || _CSPE > _cycleNum){
            _resState = "===AnalysisResult: Failed- CSPE out of range";
            return -1;
        }
        if(func->a < _minRange){
            _resState = "===AnalysisResult: Failed- func->a < minRange";
            return -1;
        }
        if(func->b > 9.5 && _functype>0){
            _resState = "===AnalysisResult: Failed- func->b > 9.5";
            return -1;
        }
    }


    // 手动计算CSDM
    double sdm = -1, CSDM;
    for(double c = _CSPE; c < _CFDM; c += 1){
        if(func->diff2(c) > sdm){
            sdm = func->diff2(c);
            CSDM = c;
        }
    }_CSDM = CSDM;
    for(double c = CSDM-0.6; c < _CSDM + 0.6; c += 0.05){
        if(func->diff2(c) > sdm){
            sdm = func->diff2(c);
            _CSDM = c;
        }
    }


    _threshold = func->f(_CSPE);

    // 扩增效率
    _E = calEfficiency();
    _slope = - 1 / log10(_E);

    return 0;
}



#pragma endregion CALCULATION}
